First-order phase transition and anomalous hysteresis of Bose gases in optical lattices 
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We study the first-order quantum phase transitions of Bose gases in optical lattices. A special 
emphasis is placed on an anomalous hysteresis behavior, in which the phase transition occurs in a 
unidirectional way and a hysteresis loop does not form. We first revisit the hardcore Bose-Hubbard 
model with dipole-dipole interactions on a triangular lattice to analyze accurately the ground-state 
phase diagram and the hysteresis using the cluster mean-field theory combined with cluster-size 
scaling. Details of the anomalous hysteresis are presented. We next consider the two-component 
and spin-1 Bose-Hubbard models on a hypercubic lattice and show that the anomalous hysteresis 
can emerge in these systems as well. In particular, for the former model, we discuss the experimental 
feasibility of the first-order transitions and the associated hysteresis. We also explain an underlying 
mechanism of the anomalous hysteresis by means of the Ginzburg-Landau theory. From the given 
cases, we conclude that the anomalous hysteresis is a ubiquitous phenomenon of systems with a 
phase region of lobe shape that is surrounded by the first-order boundary. 

PACS numbers: 03.75.Lm, 03.75.Mn, 05.30.Jp 



I. INTRODUCTION 

Systems of ultracold gases confined in optical lattices 
have been often used as quantum simulators for the stud- 
ies of strongly interacting many-body physics [l(. It is 
of great advantage that many essential features, includ- 
ing the statistics of particles, the optical-lattice depth, 
the filling factors, and the interparticle interactions, are 
widely controllable. The ultimate goal of the quantum 
simulation is to emulate the complex models that can- 
not be accurately solved with currently available theoret- 
ical approaches, such as the Fermi- Hubbard model and 
the quantum spin models on a kagome lattice. However, 
to that end, it is also important to examine in advance 
the performance of the optical-lattice systems as simula- 
tors by comparing experiments with theories in solvable 
models. Trotzky et al. has indeed presented quantita- 
tive comparison between the experiment with spinlcss 
Bose gases in a cubic optical lattice and the quantum 
Monte Carlo (QMC) analyses on the corresponding Bose- 
Hubbard model to show that they agree well with each 
other regarding phase transitions from supcrfluid (SF) 
to normal at finite temperatures upon approaching the 
Mott transition point Q. 

The recent development of experimental techniques 
has offered possibilities to simulate more complicated and 
richer systems. For example, mixing different hyperfine 
[a-Hlj. atomic species 
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states [3H11IJ. atomic species [12H14I] . or isotopes 
optical lattices has added pieces to access a variety of 
quantum phases and phase transitions. Moreover, a long- 
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range and anisotropic interaction between particles has 
been introduced by the creation of degenerate gases of 
atoms with stron g m agnetic dipole-dipole interactions, 
such as chromium [la ] , dysprosium [171 ] , and erbium [181 ] . 
Experimentalists have also attempted to realize the quan- 
tum degeneracy in gases of heteronuclcar polar molecules 
with stron ger dipolar interactions, such as KRb Qjl, [2CJ 
and LiCs [21j . It is also worth noting that even the 
lattice geometry is flexibly controllable; two-leg ged lad- 
der @, i, [H-pl , triangular (HHl, hone y comb l^MM , 
and kagome [3CJ optical lattices have been created in re- 
cent experiments in addition to the standard hypercubic 
lattices [3lT 33j . Such rapid expansion of the range of 
application demands further sophistication of the perfor- 
mance of quantum simulators. 

More specifically, it has been predicted that there ex- 
ist many quantum phase transitions of first order in 
some of the systems described above [34Ho7| . Hence, the 
quantum simulator has to be able to distinguish a first- 
order (discontinuous) phase transition from a second- 
order (continuous) one in order to map out correctly the 
phase diagrams of those systems. Furthermore, since the 
first-order transition is one of the fundamental subjects in 
thermodynamics and statistical physics, the observabil- 
ity of the accompanying phenomena such as hysteresis 
can be a touchstone to test the performance of a quan- 
tum simulator based on trapped atomic/molecular gases 
in optical lattices. 

The first-order quantum phase transitions in the Bose- 
Hubbard systems are interesting in part because they 
occasionally exhibit an anomalous hysteresis behavior as 
predicted for the melting transition of the hardcore Bose- 
Hubbard model with long-range interactions on a trian- 
gular lattice (521 ] . In this anomalous hysteresis, the phase 



transition occurs only unidircctionally; while the solid or- 
der melts when the chemical potential varies, the system 
in the SF state cannot be solidified again in the inverse 
process. However, comprehensive understanding of this 
phenomenon has not been achieved yet in the sense that 
the following two questions are open. First, what are 
the necessary conditions for the anomalous hysteresis to 
emerge? This question can be inductively addressed if 
one shows a few other examples. Second, can the hystere- 
sis be described within the Ginzburg-Landau theory for 
the first-order transitions? Answering the second ques- 
tion is important because the anomalous hysteresis seem- 
ingly contradicts with the Ginzburg-Landau theory that 
always predicts the spinodal points upon varying mono- 
tonically the coefficient relevant to the transition. 

In this paper, we study quantum phases and phase 
transitions of Bose gases in optical lattices, with a cen- 
tral focus on the first-order transitions. We consider the 
systems of (i) dipolar hardcore bosons on a triangular 
lattice and (ii) multi-component bosons on a hypcrcubic 
lattice. Both systems are known to exhibit a first-order 
phase transition. We discuss the first-order transition 
phenomena of these systems, including the possibility of 
experimental observations and the expected characteris- 
tic hysteresis behavior, in a comprehensive manner. 

In the previous work on system (i), Yamamoto et 
al. mainly analyzed a simple model in which the infinite- 
range dipolc-dipole interaction is truncated to be only 
nearest- neighbor (NN) interaction [52|]. We first recon- 
sider in more detail the first-order transitions for the 
model with full dipolc-dipole interaction using a larger- 
size CMF method (with up to 15-site clusters) and the 
cluster-size scaling. We find that the untruncated model 
can also exhibit the anomalous hysteresis. Besides, we 
quantitatively evaluate the phase boundaries through the 
size-scaling of the CMF data and predict the presence of 
several additional phases that have not been reported in 
the previous QMC simulations |51| . 

As for system (ii), we specifically investigate the first- 
order phase transition from SF to the Mott insulator 
(MI) in the two-component and spin-1 Bosc-Hubbard 
models within the site-decoupling mean-field approxima- 
tion. We show that this first-order transition exhibits 
the anomalous hysteresis behavior as well as in system 
(i). The two-component model has been already emu- 
lated in several previous experiments 0-0, [l2| • Having 
these experiments in mind, we discuss what conditions 
are required for observing the hysteresis associated with 
the first-order transition. Especially, we calculate the 
transition temperature from SF to normal phase upon 
approaching the Mott transition, in order to reveal the 
requirements regarding the temperature and the lattice 
depth. Moreover, taking advantage of the simplicity of 
the two-component model with respect to the first-order 
transition, we construct the Ginzburg-Landau theory for 
the anomalous hysteresis to elucidate its physical mech- 
anisms. Extracting common features from the obtained 
examples, we finally argue that the anomalous hysteresis 



emerges universally when a phase region of lobe shape is 
surrounded by a first-order boundary. 

The remainder of the paper is organized as follows. 
In Sec. HIl we investigate the hardcore Bosc-Hubbard 
model with long-range dipole-dipole interaction. The 
ground-state phase diagram obtained from the cluster- 
size scaling of the CMF data is compared with the pre- 
vious QMC data in Sec. Ill Al and the hysteresis behav- 
ior accompanying first-order transitions is discussed in 
Sec. lIIBl In Scc. lIIII we analyze the two-component Bosc- 
Hubbard model. We discuss the first-order SF-to-MI 
transitions induced by varying the hopping amplitude or 
the chemical potential. After the site-decoupling mean- 
field theory is introduced in Sec. IIII Al the cases of equal 
and unequal intra-component interactions are studied in 
Sec. IfflBl and Sec. UrTCl In Sec. IBID! we construct the 
Gintzburg-Landau theory for the anomalous hysteresis 
phenomenon. In Scc. lIVI the spin-1 Bose-Hubbard model 
is analyzed. The results are summarized in Sec. [Vj 



II. DIPOLAR BOSONS IN A TRIANGULAR 
LATTICE 



We consider a dipolar Bose gas loaded into a deep tri- 
angular optical lattice and describe the system with the 
following Bosc-Hubbard model, 

H = -J^2(a]ai + H.c.) + -^Vjifijhi - fx^iij, (1) 
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where a! and hj = alcij are the creation and number 
operators of the hardcore bosons at site j, J denotes the 
hopping amplitude between NN pairs, and fi the chemical 
potential. Here we take the hardcore boson limit, which 
means that two or more bosons are not allowed to occupy 
the same site due to the strong on-site interaction. We 
assume that the dipolc moments are polarized by the 
external field in the direction perpendicular to the lattice 
plane such that the dipole-dipole interaction Vji is well 
approximated as [58| 
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where d is the lattice spacing and Tj = (j x d,j y d) with 
integers j x and j y is a lattice vector at site j. 

Since the dipolc-dipole interaction rapidly falls off as 
the inverse cube of the distance \rj — rj|, several essen- 
tial properties may be captured using a simplified model 
with keeping only the NN repulsion V. Such a truncated 
model has been the subject of intensive study [43450l l52- 
|54| due to its simplicity. The outline of the ground-state 
phase diagram for this model has been given by the ear- 
lier MF study HH and QMC calculations [IH . Recently, 
analyses with a large-size CMF approach have pointed 
out that the SF-to-SS transition should be of first or- 
der [52| , and this has been reconfirmed by the latest QMC 



analyses [53|, l54[. Thus, the complete ground-state phase 
diagram is well established for the truncated model. 

For the case of full dipole-dipole interaction given in 
Eq. ((2J, only a few previous works [5l|,l52f have calculated 
the ground-state phase diagram. The QMC analyses in 
Ref. [5l| have shown that there exist the two solid phases 
with one-third and two-third fillings, and the supersolid 
phase in between the two solids, as in the case of the 
truncated model. Neither additional solids nor supcr- 
solids were found. However, the error bars attached to 
the phase boundaries are rather large in the QMC data 
of Ref. [5l| so that some of important properties of the 
phase diagrams might have been missed. Hence, more 
careful analyses are needed for further refinement of the 
phase diagram. 



A. Ground-state phase diagram 

We re-examine the ground-state phase diagram of the 
dipolar Bose-Hubbard model by means of the large-size 
CMF method and the cluster-size scaling analysis [571 ]. 
We use the series of the clusters that consist of Nq =3,6, 
10, and 15 sites. The shape of each cluster is depicted in 
TablcQ] In the CMF method, we perform exact diagonal- 
izations of the cluster system with mean-field boundary 
condition. The values of the mean fields are determined 
by solving the CMF self-consistent equations (see Ref. [57] 
for more details). Furthermore, we carry out a cluster- 
size scaling of the CMF data with the scaling parameter 
A defined by Nb/{Nc x z/2), where JVb is the number of 
bonds within the cluster and z = 6 is the coordination 
number of the triangular lattice. 

In the CMF calculation, one assumes background sub- 
lattice structures that are expected to emerge from the 
given lattice geometry and the nature of interactions |57| . 
In the simple truncated model, only the states with 
\/3 x \/3 sublattice structure [Fig. QJa)] emerge as a 



TABLE I: The series of the clusters used in the CMF calcu- 
lations. The values of Nb and A are also listed. 
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FIG. 1: (a) A three-sublattice v3 x \/H structure and (b) 
two types of four-sublattice structures assumed in the present 
calculations. The lattice sites with the same number belong 
to the same sublattice. 



natural result of the interplay between the NN density- 
density repulsion and the geometry of the triangular lat- 
tice [44l - [50l I52M54I ] . In the case of full dipole-dipole inter- 
action, more complicated patterns that have a larger unit 
cell are also expected to be formed since the interaction 
can act on longer-distant boson pairs [59j . Therefore, 
we take into account two types of four-sublattice pat- 
terns shown in Fig. QJb) m addition to the basic three- 
sublattice -\/3 x y/3 structure. 

Figure [2] shows the ground-state diagram of the dipo- 
lar Bose-Hubbard model given in Eq. (fTJ) obtained by the 
CMF method with the 15-site cluster (CMF-15). The 
phase diagram is symmetric around p/V = (p/V)o = 
J2i Vji/2V « 5.517, reflecting the particle-hole symme- 
try of the hardcore boson system. We see the solids with 
the filling factor p = 1/4, 1/3, 1/2, 2/3, and 3/4, and 
the two types of SS, which we name SSI and SS2, as 
well as the standard SF state. While the solid phases 
form a superlattice pattern in the local density {fij), i.e., 
a crystalline order, the SF state is identified by the or- 
der parameter \& = Ei(%)A^> where M denotes the 
number of lattice sites. The SS phases simultaneously 
possess the SF order and the crystalline order. To de- 
termine the boundary lines of first-order transitions, we 
used the Maxwell construction in (J/V,x)-plane, where 
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Solving the CMF self-consistent equations with the 
ansatz of sublattice structures shown in Figs. [TJa) 
and []Jb) , we obtain the crystalline order of each solid 
or SS state depicted in the lower panels of Fig. [2] As 
shown in the phase diagram, the solids with p = 1/3 and 
2/3, and SSI are present at relatively large J/V, as in 
the truncated model |44l. l45l. |52H54| . However, there also 
emerge the solids with p = 1/4, 1/2, and 3/4, and SS2 
due to the long-range nature of the dipole-dipole interac- 
tion. Moreover, many other solid and SS phases that 
cannot be described by the three-sublattice and four- 
sublattice ansatz should be found for smaller values of 
J/V, as in the case of a square lattice |59j . 

In the previous work (52|, we presented the single-site 
mean-field (MF) and CMF-10 results of the ground-state 
phase diagram (see Fig. 4 in Ref. 1 521 ). Compared to 
them, the solid and SS regions in Fig. [5] arc shifted to the 
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FIG. 2: Ground-state phase diagram of dipolar bosons in a 
triangular lattice obtained by the CMF-15. Second- and first- 
order phase transitions are indicated by the thin and thick 
lines, respectively. Lower panels: Schematic pictures of the 
density patterns in each phase. Because of the particle-hole 
symmetry, we show only the phases on the low-density side 
of the phase diagram. 



side of small J/V and especially the region of SS2 phase 
is much narrower. This means that the results of the 
phase boundaries are not completely converged. In order 
to include systematically the effects of quantum fluctu- 
ations, we carry out the cluster-size scaling of the CMF 
data with different sizes of clusters. We perform linear 
fits of the data obtained from the three largest clusters, 
namely Nq = 6, 10, and 15. The advantage of this scal- 
ing method is that usually the series of CMF data are 
well fitted to the linear function even if the clusters are 
not quite large [57]]. This is attributed to the fact that 
the CMF method treats infinite-size systems by setting 
the mean-field boundary condition despite using a finite- 
size cluster. This point is an important difference from 
the standard exact diagonalization with open or periodic 
boundary condition. Indeed, this scaled CMF procedure 
has successfully generated a quantitatively reliable result 
for the phase boundaries of the dipolar model on a square 
lattice [571 ]. 

Figure [3] shows the resulting scaled data of the phase 
boundaries in the ( J/V,ju/V)-plane. Although the solid 
and SS regions are shifted further to the side of small 
J/V, the relative location of the phases does not change 
even in the limit of Ac — > oo except for the following 
point: The calculations with up to Nq = 15 clusters show 
that a narrow but finite region of the SS2 phase exists 
around half filling and the transition from the low-filling 
SSI to high-filling SSI state occurs continuously through 
the intermediate SS2 state (see Fig. [5]). However, the SS2 
region rapidly shrinks with increasing the cluster size and 




0.06 0.08 0.10 0.12 0.14 

JIV 

FIG. 3: (color online) The result of the cluster-size scaling for 
the phase boundaries. For comparison, we al so p lot the QMC 
data of Ref. [5l| (green dashed lines; see Ref. [5l| for details). 



completely disappears in the limit of Nq — > oo. Hence, 
the low-filling SSI state is directly transformed into the 
high-filling SSI state with a finite jump in, e.g., the filling 
factor at the particle-hole symmetric line fx/V = (p,/V)o. 
This first-order nature of the SS1-SS1 transition has been 
predicted also for the truncated model {HEtT. |49H52J| . 

As shown in Fig. Ufa), the linear fits of the CMF 
data for the SF-SS1 and p = 1/3 solid-SSI boundaries 
are fairly good. In fact, we can see in Fig . [3] that the 
scaled values and the QMC data of Ref. [5l| are compati- 
ble within the error bars. However, in contrast to Ref. l5ll 
the SSI phases are found not only in between the p = 1/3 
and 2/3 solids but also for n < 1/3 and n > 2/3. In ad- 
dition, our scaled CMF result predicts the existence of 
solid states with p = 1/4, 1/2, and 3/4 in the plot range 
of the figure, although these states have not been de- 
tected in the QMC analysis for J/V > 0.045 [p. As for 
the p = 1/2 solid, the scaled phase boundary is reliable 
because the data of the three largest clusters almost lie 
in a straight line as seen in Fig. 01b). However, the lin- 
ear fitting may not make a very good estimation for the 
p = 1/4 solid boundary. Therefore, more careful and pre- 
cise analyses are demanded to corroborate the presence 
of the new phases. Such analyses may be able to be done 
with the use of the latest QMC algorithm developed in 
Ref. [60| for treating efficiently long-range interactions. 



B. First-order phase transitions and anomalous 
hysteresis 

The system of dipolar bosons in a triangular lattice ex- 
hibits first-order phase transitions indicated by the thick 
lines in Figs. [5] In Fig. [SJ the filling factor p shows finite 
jumps at the first-order transition points. We focus on 
the quantum melting transition from solid to SF around 




FIG. 4: (a) Examples of the cluster-size scaling of the CMF 
data for the phase boundaries J c /V between the SF and SSI 
phases (open symbols), and between the SF and p — 1/3 solid 
phases (closed symbols), (b) Cluster-size scaling for the tips 
of the p — 1/4 and 1/2 solid lobes. 



the plateau with p = 1/3. In Ref. [52J, it has been shown 
for the truncated model that the melting transition shows 
an anomalous hysteresis behavior. We establish here that 
the anomalous hysteresis can occur when the full dipole- 
dipole interaction is taken into account. 

In Fig. Eta), we plot the solid order parameter of the 
\/3 x \/3 sublattice structure, 
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(3) 



as a function of p/V for J/V — 0.14. The thick solid lines 
represent the value of pq at the ground state. Points a 
(a') and b (b') correspond to the first-order transition 
(equal-energy) points. The curve of pq at the ground 
state exhibits finite jumps at these points as well as the 
curve of p in Fig. \5\ The CMF method can calculate also 
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FIG. 5: The filling factor p as a function of p/V for J/V = 
0.14 (obtained by CMF-15). The first-order transition points 
are located at the discontinuous jumps in p. 
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FIG. 6: (a) The solution curves of CMF-15 for the solid or- 
der parameter |pq| as a function of p/V at zero temperature 
and J/V = 0.14. The thick solid, thin solid, and dashed lines 
represent the ground, metastable, and unstable states, respec- 
tively, (b) and (c) The transition processes in the anomalous 
hysteresis. 



metastable and unstable solutions. The solution curve 
of pq for metastable (unstable) states is shown by the 
thin solid (dashed) lines. The total solution curves form 
one line of the SF solution (e-/) and one closed loop con- 
sisting of the solid and SSI solutions (a-b-y-x-a) . Such 
topology of the solution curves is in stark contrast to 
conventional first-order phase transitions, e.g., a typical 
liquid-solid transition. The formation of the separated 
solution curves originates from the re-entrance of the 
SF phase near the tip of p = 1/3 (or 2/3) solid lobes 
in the phase diagram upon sweeping the value of p/V. 
Of importance is that the ground-state SF solutions in 
the small and large p/V sides are connected through the 
metastable solutions, which means that the SF state re- 
mains (at least locally) stable over the entire range of 
p/V in Fig. IH^a). This fact leads to a characteristic hys- 
teresis behavior as described below. 

Let us assume that we initially prepare a solid state 



which is located at point i in Fig. HJa) as an initial state. 
When fi/V decreases, the system undergoes a melting 
transition to the SF phase though the transition path 
i — > x — > x' — > e [upper panel, Fig. GD^ 3 )]. The solid 
state can remain metastable even after fi/V exceeds the 
first-order transition point b, and the melting of solid or- 
der occurs at the metastability limit of the SSI phase 
(point x). Next, let us consider the reverse process (so- 
lidification) from the SF state at point e. In this case, 
the system remains in the SF phase upon increasing fi/V 
through e — > / [lower panel, Fig. E][b)]. This is indeed 
attributable to the absence of the metastability limit of 
the SF phase. In a usual hysteresis, the reverse transition 
also can occur and the transition path forms a hysteresis 
loop. However, in this anomalous hysteresis, the tran- 
sition can occur only unidirectionally from the solid to 
SF state, and the SF state cannot be solidified again by 
varying fi/V. A similar anomalous hysteresis behavior 
also occurs when fi/V first increases and then decreases 
[see Fig. GH^c)]. 



It is practically difficult to control the global chemi- 
cal potential directly by external manipulation on cold 
atomic gases. However, one could vary the local chem- 
ical potential in the following way. In actual experi- 
ments of ultracold gases, there exists a trap potential, 
e.g., Vt(r) — mu) 2 \r\ 2 /2 in addition to the optical lattice 
potential. Within the local-density approximation, the 
local chemical potential is written as p,j = fi — Vt (tj). 
Therefore, we suggest that the anomalous hysteretic be- 
havior could be realized by manipulating the frequency 
ui to control flj at the trap center. 



Figures Eta) andjTJb) plot the metastability- limit line 
of the SF phase in the phase diagram obtained by CMF- 
15 and the cluster-size scaling. The anomalous hystere- 
sis can be observed when the value of J/V lies within 
the range sandwiched between the two vertical dotted 
lines, in which the SF state is always stable for any value 
of fi/V. In this case, the lobe region consisting of the 
normal solid and SSI phases is surrounded by the first- 
order transition boundary with the SF phase. In the 
presence of such a lobe region with the first-order transi- 
tion boundary, the metastability limit of the surrounding 
phase (the SF phase in this case) is necessarily located 
inside the lobe as in Fig. and therefore there must be a 
finite region where the surrounding phase is always stable 
for any value of the parameter in the vertical axis {fi/V 
in this case). This implies that the anomalous hysteresis 
is not a unique phenomenon in this dipolar system but a 
common feature of the systems that have the geometry 
of the phase diagram mentioned above. In the following 
sections, in order to make this implication more convinc- 
ing, we will provide a few other examples of the systems 
that exhibit the anomalous hysteresis. 
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FIG. 7: The enlarged views of the phase diagrams obtained by 
the CMF-15 calculation (Fig. [2]) and the cluster-size scaling 
[Fig.[3jb)]. The dashed curves indicate the metastability limit 
of the SF phase. On the right side of the curve, the SF state 
is always (meta)stable. 



III. BOSE-BOSE MIXTURES IN A 
HYPERCUBIC LATTICE 

In this section, we consider the system described by 
the following two-component Bose-Hubbard model [61| 
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where b- a creates a boson of the particle type a 
at site j of a D-dimensional hypercubic lattice, fij >a = 
b] a bj :a , t a is the NN hopping amplitude, U a the on- 
site intra-component interaction, XJ\i the on-site inter- 
component interaction, and fi a the chemical potential. 
This model system could be realized using a mixture of 
two different types of bosons loaded into a sufficiently 
deep optical lattice 0-0 13 ■ The strength of U±2 can 
be controlled experimentally by using Feshbach reso- 
nances I5j . |62l . 1(33 1 or component-dependent optical lat- 
tices |7l.l64|. 

Previous theoretical works have extensively investi- 
gated the ground-state phase diagrams of Eq. (TO for dif- 
ferent parameter regions and dimensionality [3414371 l65l — 
|74| . In particular, it was found that the transition from 



SF to MI with even fillings can be of first order [35l — 
|37| . Since the MI region takes the lobe shape mentioned 
in the previous section, the anomalous hysteresis is ex- 
pected to occur. In the following we will show that this 
is indeed the case for both equal and slightly-unequal 
intra-component interactions. For observing the SF-to- 
MI transition of first order, the two-component system 
is advantageous over the other systems in the sense that 
it has been already realized in several experiments [j- 
0jLL2J]- To stimulate further experimental efforts, we esti- 
mate the temperature and the lattice depth at which the 
first-order transitions can be clearly observed. Finally 
in Sec. IIII D[ we present the Ginzburg-Landau theory for 
the first-order SF-to-MI transitions to explain the anoma- 
lous hysteresis within the Ginzburg-Landau framework. 



A. Site-decoupling mean-field approximation 



We analyze the Hamiltonian in Ec 
(single-site) MF approximation [75M7 
hopping term as 



using the simple 
Decoupling the 



b]Ji,a « (b}Jb ha + (b La )b} a -(b] a )(b^), (5) 

we approximate the system by the sum of M identical 
MF Hamiltonians: H w £\ Hf F with 



H 



MF 



o=l,2 



(0 



J j,a V3,cc 



flaUj 



U-lihjAUjo- 



(6) 



Here the sum X)m runs over z — 2D nearest- neighbor 
sites of site j and M is the number of lattice sites. The 
mean field cf>j ta = (bj ta ) = (bj a ) plays the role of the 
SF order parameter. Here, we take the order parameters 
cj)j t i and cj)j t 2 to be real without loss of generality. Under 
the assumption of the spatial homogeneity (f>j t0l = <f> a , the 
many-body lattice problem is reduced to a set of indepen- 
dent single-site problems with the effective Hamiltonian 
7Jj vIF and therefore we drop the site-labeling subscripts 
as 

rVMF . . r>MF 



= ^2 \-zt Q (f> Q ( b a 



+yi(n« - 1) 



l^a^c 



Ul2nifl2. 



(7) 



The value of the SF order parameter 
that the free energy per site 

. , 1. 



is determined so 



F( 



P' 



:Tr 



( exp [-M MF : 



(8) 



takes a minimum value with respect to </>i and 02- Here 
f3 = l/(k&T) is the inverse temperature. Around a first- 
order phase transition, the energy function has two or 
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FIG. 8: The ground-state phase diagram of the isotropic two- 
component Bose-Hubbard model within the MF approxima- 
tion for U12/U — 0.9. Second- and first-order phase transi- 
tions are indicated by the thin and thick lines, respectively. 
The phase boundaries are identical to those in Fig. 1(a) 
of Ref. [33|. The dashed (dash-dotted) lines represent the 
metastability limits of SF (MI) phase, and the dots mark the 
tricritical points, where the transition changes from first to 
second order. The inset is an enlarged view of the region 
around the tip of the p = 2 MI lobe. 



more local minima, one of which corresponds to the glob- 
ally stable state and the others are mctastable states. 
Additionally, F(<j>i,<J)2) has several maxima or saddle 
points corresponding to unstable stationary states. In- 
stead of the minimization of the free energy, we can also 
obtain <p\ and <p2 by directly solving the set of two self- 
consistent equations 

4> a = Tr (b a cxp[-/3H MF ]] /Tr f exp[-/3H MF ]) (9) 



with a — 1 and 2. 

It is known that the site-decoupling MF method can 
produce a qualitatively correct result for phase transi- 
tions between SF and MI phases [75|, UM ■ This method 
becomes exact in the limit of infinite spatial dimen- 
sions. Therefore, the MF Hamiltonian H MF is expected 
to give a good approximation to binary Bose mixtures 
in a three-dimensional (simple-cubic) optical lattice with 
z = 6. Another essential advantage for this study is that 
metastability analyses can be easily performed on the 
basis of the energy function F(<j>i,<j>2). 



B. Equal inter-component interactions 

In this section, we focus on the simple case of t\ = t% = 
t, U\ = U2 = U > 0, and \i\ = n% = /1. Because of the 
1 •(-> 2 exchange symmetry of the Hamiltonian, we take 
0i = <f>2 — <t>- We briefly review the MF ground-state 
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FIG. 9: (a) The MF solution curve (left) and the expected 
hysteresis loop (right) in the SF-to-MI transition for p/U = 
1.314 at zero temperature. In the left panel, the thick solid, 
thin solid, and dashed lines represent the ground, metastable, 
and unstable states, respectively. The vertical dotted line 
indicates the equal-energy point (the first-order transition 
point), (b) The U12/U dependence of <5h and <5 a - 



phase diagram for U12/U = 0.9 shown in Fig. [5J where 
the phase boundaries have been obtained in Fig. 1(a) of 
Ref. [37| . First-order phase transitions arc found near 
the tips of the MI lobes (thick lines) unlike the SF-to-MI 
transition in the single-component Bosc-Hubbard model. 
We discuss the range of the hysteresis region and the 
temperature effects in the first-order transition between 
the SF and p = 2 MI phases because they are essential for 
the visibility in actual experiments. We consider the SF- 
to-MI transition upon cycling (increasing and decreas- 
ing) zt/U. The corresponding experiment has been done 
for single-component Bose gases by tuning the depth of 
the lattice potential 0, l3ll - l33| . In the case of Bosc-Bosc 
mixtures, the system is expected to show a hysteresis be- 
havior in the SF-to-MI first-order transition. Figure [§ta) 
shows the solution curve of Eq. © and the expected 
hysteresis accompanying the phase transition between 
SF and MI p = 2 for p/U = 1.314. In this case, the 
transition process forms a conventional hysteresis loop 
(right panel). Unlike the anomalous hysteresis discussed 
in Sec. UH the transition occurs bidirectionally between 
the SF and MI phases. The width 6^ indicated by the 
arrow in the inset of Fig. [3] corresponds to the maximum 
size of the emergence region of the hysteresis behavior. 
We see that the hysteresis region is sufficiently sizable 
(w 13% of the size of the p = 2 MI lobe in zt/U). As 
shown in Fig. HJb), the width Sh gets wider as U12/U is 
increased. If the inter-component repulsion is not suffi- 
ciently large (U12/U < 0.68), the hysteresis region van- 
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FIG. 10: Finite-temperature phase diagrams of the two- 
component Bose-Hubbard model for (a) U12/U = and (b) 
U12/U = 0.9. The chemical potential p is tuned to be the 
value at the tip of the p — 2 MI lobe; p — 0.414 for (a) and 
p = 1.314 for (b). Second- and first-order phase transitions 
are indicated by the thin and thick lines, respectively. In ad- 
dition to the SF transition temperature T c /U, we plot the 
value of the Mott gap Ami/U as a function of zt/U. Lower 
panel of (b): the enlarged view of the upper panel with adding 
the metastability limits of the SF (dashed line) and MI (dash- 
dotted line) phases. 



ishes and the SF-to-MI transition becomes second order 
at the entire boundary of the p = 2 MI lobe [37| . 

Let us discuss the effects of finite temperatures on the 
first-order nature of the SF-to-MI transition. Figure [TOl 
shows the SF transition temperature k-oT c /U and the 
Mott gap Ami/U as functions of zt/U, where U12/U = 
(a) and 0.9 (b). When U12 = 0, the Hamiltonian is com- 
pletely separated into two independent single-component 
Bose systems. Therefore, the result shown in Fig. [TUTa) 
is identical to the SF-to-MI transition at unit filling in 
the single-component Bose-Hubbard model. At zero tem- 
perature, the system undergoes a continuous quantum 
phase transition from the SF to the MI phase at the crit- 
ical ratio of the hopping amplitude to the interaction 
strength. The MF approximation leads to the critical 
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FIG. 11: (a) The solution curves for the SF order parameter <j> 
as a function of \ijTJ at zero temperature and zt/U — 0.153. 
The thick solid, thin solid, and dashed lines represent the 
ground, metastable, and unstable states, respectively, (b) 
and (c) The transition processes in the anomalous hysteresis. 



value (zt/U) c = 1/5.83 (ZZllZa]- The transition temper- 
ature T c and the Mott gap Ami simultaneously vanish 
at the quantum critical point (QCP) [2|. Notice that 
the Mott gap Ami corresponds to the width of the MI 
lobe in (i at zero temperature. In a strict sense, the 
MI phase exists only at T = 0. However, for low tem- 
peratures, Mott-like features still remain in the normal 
fluids, where the filling factor deviates only slightly from 
integer values. When U12/U > 0.68, the situation dras- 
tically differs from the single-component (or U12/U = 0) 
case. As shown in Fig. [TWb). the SF-to-MI transition 
for low temperatures becomes of first-order and thus the 
QCP vanishes. In the region of the first-order transi- 
tion, Ami is given by the width between the upper and 
lower metastability limits of the MI phase in the phase 
diagram (see Fig. [8]) . The first-order nature of the SF- 
to-MI (-normal) transition remains until T ps 0.08U for 
U12/U = 0.9. 

We next consider to vary li/U at fixed zt/U such that 
the anomalous hysteresis occurs. We see in Fig. [8] that 
the system undergoes a re-entrant first-order transition 
from SF to MI, and back to SF near the tip of the Mott 
lobe with even fillings. For a certain value of zt/U, the 
SF phase remains stable for any fi/U, and the hystere- 
sis process shows the anomalous unidirectional behavior. 
Figure QTJa) shows the solution curves of Eq. @ as a 
function of li/U around the transition between the SF 
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FIG. 12: The curve of the MF energy measured from the en- 
ergy of the stable SF state in the re-entrant first-order tran- 
sition between SF and MI. The thick solid, thin solid, and 
dashed lines are the global minimum, local minimum, and 
maximum values of the energy landscape, respectively. Each 
point corresponds to the point with the same letter in Fig. 1111 



and p = 2 MI phases for zt/U = 0.153 at zero temper- 
ature. We see that the solution curves are completely 
separated into one closed loop and one line unlike a con- 
ventional first-order transition [compare with Fig. IHta)] . 
Thus, the hysteresis exhibits a unidirectional behavior 
upon cycling li/U . 

The transition process is similar to the one of the 
triangular-lattice dipolar system discussed in Sec. HU Let 
us start from an initial MI state at point i. If li/U is 
decreased and then increased, the transition from MI to 
SF can occur (i — > x — > x' — > e) but the system does not 
return to MI (e —> /) as shown in Fig. UTT b). The unidi- 
rectional behavior is also obtained upon increasing first 
and then decreasing \x/U [Fig. Illf c)]. The energy curve 
in the anomalous hysteresis has a characteristic shape 
where a line pierces a closed loop as shown in Fig. [T2] 
This is obviously different from the one in the conven- 
tional first-order transition, which forms a so-called swal- 
lowtail structure [80|. The anomalous behavior can be 
seen when zt/U lies within the width <5 a in the inset of 
Fig. where the SF phase is stable for any ll/U . As 
shown in Fig. HJb), the width d a exhibits a maximum at 
Uu/U ra 0.98. 



C. Unequal intra-component interactions 

While we assumed U\ = U<z for theoretical simplicity in 
Sec. IIIIB| this condition does not hold in actual experi- 
mental systems of Bose-Bose mixtures. We here consider 
the case of U\ 7^ Ui to discuss the first-order SF-to-MI 
transitions in a realistic situation. Specifically, we sup- 
pose a binary mixture of 87 Rb atoms in the two hypcrfine 
states, \F = l,m,p = 1) and \F = 2,mp = —1), which 
has been created in several previous experiments [j-Q]- 
We label the former (latter) hyperfine state as a = 1 
(a = 2). Since the scattering lengths for the intra- 
component interactions are given by a\ = 100.40ob and 
<i2 = 95.00q,b [63|, |81|, U\ and U2 are slightly unequal as 
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U2/U1 = 0.9462, where Ob is the Bohr radius. We recall 
that U12 can be widely controlled using the Feshbach res- 
onance or the component-dependent optical lattice, while 
the bare scattering length is given by a\2 = 97.66ob. 
When U\ ^ U2, the SF order parameters 4>i and 4>2 can 
take different values, and it is even possible that one com- 
ponent of bosons forms a MI state while another compo- 
nent is condensed; e.g., 4>\ = and <p2 ^ [351 . We call 
this phase MIi+SF 2 at zero temperature and NFi+SF 2 
at finite temperatures, where NF means a normal fluid. 
We tune the difference of the chemical potentials /ii — /12 
so that the populations of the two components are bal- 
anced, i.e., (hi) = (77.2). The total filling factor p is con- 
trolled by the averaged chemical potential fl = ^2 a fi a /2. 

Solving the set of Eqs. (|9]) for <\>\ and <p2, we deter- 
mine the ground-state phase diagram for t\ = £2 = t, 
\J2IVx = 0.9462, and U l2 /Ui = 0.9. We sec in Fig. QJ 
that the slight difference of the intra-component interac- 
tions hardly changes the phase diagram. The MI1+SF2 
phase emerge only as a metastable state in a very small 
region near the tip of the metastability limit of each MI 
phase (see the inset). Besides, Fig. [i"4l shows that the 
finite-temperature phase diagram also remains almost 
the same as the case of U\ = U2 in Fig. fTUTb). The 
only one noticeable change is the emergence of a narrow 
region of the NF1+SF2 phase at finite temperatures. 

Since there remain the first-order SF-to-MI transitions 
and the hysteresis, it is expected that one can simulate 
them by using the binary mixture of 87 Rb atoms in op- 
tical lattices 0-0] • As seen in Fig. [TU the temperature 
range where the first-order nature of the SF-to-MI transi- 
tion could be clearly observed is given by k-gT < 0.026^ . 
The hysteresis loop at a fixed chemical potential is ex- 
pected to be formed in the range of 0.14 < zt/Ui < 0.163 
(<5h in the inset of Fig. [j"3"|) . The anomalous hysteresis be- 
havior upon cycling the value of n/Ui is expected to oc- 
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FIG. 14: The same as in Fig. [lOjb) for U2/U1 = 0.9462 and 
U12/U1 = 0.9. The chemical potential fi is fixed to be 1.305. 
We plot the Mott gap Ami/U for the a — 1 bosons in addition 
to the finite-temperature phase boundaries and metastability 
limits. 



cur when the ratio zt/U lies within 0.14 < zt/Ui < 0.145 
(<5 a ). In the experiment with an optical lattice, the ratio 
t/U\ can be tuned by manipulating the maximum poten- 
tial depth Vq. To connect the experimental parameters 
to the Bose- Hubbard parameters, we use the following 
formulae 18211 . 
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FIG. 13: The same as in Fig. [8] for U2/U1 = 0.9462 and 
Um/Ui = 0.9. 



U 



kaE r 



Vb 

E r 



3/4 



(11) 



where E r is the recoil energy. Assuming a simple-cubic 
optical lattice with a lattice constant d = 7r/fc = 532 nm 
and the scattering length a = a\ = 100.40<2b, we obtain 
the required tuning ranges for observing the hysteresis 
loop as 13.7 < V /E r < 14.3. This level of control- 
lability of Vo/E r can be achieved by the current tech- 
niques Q, in which Vo/E r can be tuned at least in 0.25 
increments. To observe the anomalous hysteresis, the 
lattice depth has to be tuned to be within the range of 
14.2 < Vo/E r < 14.3, which is a challenge for future ex- 
perimental efforts. Besides the system needs to be cooled 
down to T < 0.02£/i w 0.70 nK for a clear observation of 
the hysteresis behavior. It is expected to be achievable 
in the near future, given the currently accessible temper- 
ature (w 1 nK |). 
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D. The Ginzburg-Landau description of the 
anomalous hysteresis 

In the previous sections, we have seen the anomalous 
hysteresis behavior in the two systems. Theoretically, the 
two-component Bosc-Hubbard system is much simpler to 
handle compared to the dipolar model in Sec. Ull because 
the anomalous hysteresis can be described only at the 
single-site MF level with one order parameter (j>. Taking 
this advantage, in this section we construct the Ginzburg- 
Landau theory for the anomalous hysteresis by taking the 
two-component Bose-Hubbard model with U\ = U 2 as an 
example. On this basis, we will clarify the reason why the 
hysteresis shows the anomalous unidirectional behavior. 

We construct the Ginzburg-Landau energy function by 
using the perturbative MF method [73, [83( in order to 
express the coefficients in the energy function as functions 
of the microscopic parameters in the Hamiltonian. Here, 
we assume t\ = £2 = t, U\ = U2 = U > 0, and \x\ = 
A'2 = fJ--, and take the SF order parameter 4>i = 4> 2 = 4> 
to be real without loss of generality. We divide the MF 
Hamiltonian H MF into diagonal and off-diagonal parts as 
H MF = H° + V + 2zt(j) 2 , where 



£°=x; 



and 



-pn a + —n a (n a - 1) 



V = -ztcf>^2(b a + bi). 



+ U l2 fnn 2 (12) 



(13) 



Assuming that (j> is small, we deal with V as a small 
perturbation. The eigenvalue of the unperturbed Hamil- 
tonian H° for the two-component Fock state \ni,n 2 ) is 
given as 



H°\n ll n 2 ) = E°(n 1 ,n 2 )\n 1 ,n 2 ), 



(14) 



where 



E°(ni,n 2 ) = -ju(ni +n 2 ) + — ni(ni -1) 

+—n 2 (n 2 -\) + \] 12 n 1 n 2 . (15) 

Now we focus on the first-order phase transitions which 
appear near the tips of MI lobes with even fillings p = 2m 
(m = 0, 1, 2, ... ). For even fillings, the ground state of 
H° is simply given by |i) = |m,m). The second-order 
correction to the energy is given by 



SE 2 = (i\V G V\i) + 2zt(f> 2 

X^ l(i|^|p)| 2 , 9 ,, 2 

P " 



(16) 



where we used the notations G = (E® — iJ ) -1 and 
A B = AJ2„ |p)(p|-B. Here J2 P runs over an the eigen- 
states of H° other than the initial state li) = \m, m). 



Ep = E°(ni,n 2 ) is the eigenvalue for the state |p) = 
|n 1 ,n 2 ). Note that the last term comes from the MF 
decoupling in Eq. ((5]). In the second-order perturba- 
tion process, we have only four intermediate states |p) = 
|m±l,m) and |m,m±l), and the matrix elements (i|V^|p) 
are easily calculated. As a result, we obtain 



SE 2 



2(zt) 2 m 



-fi + U(m,+ l) + U 12 m 



2(zt) 2 (m + l) 
fi-(U + U 12 )m 

a 2 4> 2 , 



2zt 



(17) 



where a 2 = a 2 (zt, U, Ui 2 , fj,, m) is the second-order coeffi- 
cient of the Ginzburg-Landau energy function. In order 
to discuss the first-order transition phenomena in this 
system, one needs to consider the terms up to the sixth- 
order in (j). The fourth-order and sixth-order corrections 
to the energy are given, respectively, by 



SEi 



(i\V G VoGoV G V\i) 



SE 2 (i\V Q G 2 V\i) 



«-K 



(18) 



and 



5E 6 = (i|^ G VoG V;G o yoGoV;G V|i) 

~6E 4 (i\V G 2 V\i) + (SE 2 ) 2 (i\V G 3 V\i) 
+SE 2 (i\V G 2 V GoV a G V\i) 
+6E 2 (i\V G V G 2 V G V\i) 
+8E 2 (i\V G V G V G 2 V\i) 
= a 6 6 . 



(19) 



Here, SE 2 = 8E 2 — 2zt(f> 2 . One has to consider up 
to 16 x 8 matrix elements of V between intermediate 
states in the calculation of the sixth-order coefficient 
^6 = o,g(zt, U, U\ 2 , p, 771), since many intermediate states 
emerge. However, for p = 2 (m = 1), the maximum size 
of the matrix reduces to 10 x 6 due to the lower limit of 
occupation (ni,ri2 > 0). 

Now we obtain the Ginzburg-Landau energy 



E{4>) = E® + a 2 (j) 2 + a 4 4 + a 6 6 



(20) 



whose coefficients are related to the microscopic system 
parameters in the two-component Bose-Hubbard Hamil- 
tonian through Eqs. (|17til9|) . The shape of the energy as 
a function of <p changes depending on the values of the 
coefficients. For the existence of a stable ground state, 
at least the highest-order coefficient a§ must be positive. 
If the second-order coefficient a 2 is negative, the energy 
function forms two minima (the wine-bottle or Mexican- 
hat shape when the imaginary axis of <j> is considered) 
shown in Fig. fTBTB. Therefore, the ground state is a SF 
state with a finite </>. Although the MI state with = 
is also a stationary point (dE(<f>)/d(f> — 0), it is unsta- 
ble. If a 2 becomes positive, the profile around = 
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FIG. 15: Schematic illustrations of energy landscapes. Here, 
the order parameter <f> is taken to be real. 
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FIG. 16: The phase boundary of the transition between the 
SF and p = 2 MI phases obtained by the perturbative MF 
method. We set U12/U = 0.9 as in Fig.|S] The dashed (dash- 
dotted) lines represent the metastability limit of SF (MI) 
phase, and the vertical lines with arrowheads mark the plot 
ranges of Figs. \TT\ and 1181 In the shaded region with width 
<5 a , the anomalous hysteresis appears when varying p/U. The 
tricritical points (dots) can be obtained from the condition 
a.2 = 0,4 = 0. 



changes into a convex shape, and the MI state becomes 
stable. If the transition is of second order, the condition 
a 2 = gives phase boundary from SF to MI. However, in 
first-order transitions, the energy shape can form three 
minima, corresponding to a MI state and two equiva- 
lent SF states, shown in Figs. HBJII) and [TSJIII) . There 
are two maxima which correspond to unstable SF states. 
The energy of the SF states at the minima is easily ob- 
tained by substituting the solution of dE(<fi)/d<fi — into 
Eq. ([20)) . If 04 < — 2y/a 2 a 6 , the energy of the SF state is 
smaller than that of the MI state, and the energy function 
forms a shape shown in Fig. [ToT lI). In contrast, the MI 
state has the lowest energy for 04 > — 2^/a 2 etg as shown 
in Fig. 1151111). When the value of 0,4 exceeds — v / 3a2&6, 
the stationary solutions with <f> ^ disappear and the 
energy profile has only one minimum at 4> — as shown 
in Fig. ll5f IV). In summary, the expansion up to sixth or- 
der in (j) can describe four kinds of energy shapes, which 
are separated by the conditions a 2 = 0, a\ = 4a 2 a 6 , and 
a\ = 3a2a6. These conditions correspond to the metasta- 
bility limit of MI, the first-order transition boundary (or 
the equal-energy point), and the metastability limit of 
SF, respectively. 

The coefficients a 2 , 04, and a^ are related to zt, U, U12, 
and n through the perturbative MF method. Therefore, 
the phase diagram in the (zt/U ,n/U) plane is easily ob- 
tained by the three conditions a 2 = 0, a\ = 4a 2 a 6 , and 
a\ = 3a2a,6. Figure [TH] shows the curves of the first- 
order transition boundary between the SF and p = 2 



MI phases and the metastability limits of MI and SF for 
U12/U = 0.9. Compared to the full MF result in Fig. |SJ 
the positions of the first-order transition boundary and 
the metastability limit of SF are slightly shifted, while 
the metastability limit of MI shows exact concordance. 
This is because higher-order terms are neglected in the 
perturbative expansion of the energy function. 

Let us see the changes in the shape of the energy 
function and the behaviors of the coefficients during the 
anomalous hysteresis induced by cycling the chemical po- 
tential n/U . First, we review the case of the conventional 
hysteresis loop. The value of the SF order parameter <fi 
is obtained as 



— a\ + \Jd\ - 3a 2 a 6 -a 4 - yja\ - 3a 2 a 6 



3a 6 



3a 6 



(21) 



from dE((f>)/d(f> = 0. The latter solution corresponds 
to the unstable SF states at the maxima of the energy 
function, which exist only for 04 < and a\ > 3a 2 a6, 
i.e., in the case of Figs. IToT lI) and lToT lII). We plot the 
solution curve as a function of fi/U for zt/U = 0.13 in 
Fig. ITTT a). The lower panels show the coefficients of the 
Gintzburg-Landau energy as functions of fi/U. The three 
conditions a 2 = 0, a| = 4a 2 a6, and a\ = 3a 2 a@ separate 
the plot into four regions. The numbers (I-IV) indicate 
that the energy profile in the region has the shape shown 
in Figs. [15] with the corresponding number. In this case, 
the shape of the energy function changes from (I) to (IV) 
as the value of \x/U increases. One finds three stationary 
solutions (stable MI and SF solutions and an unstable 
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FIG. 17: (a) The solution curve for the order parameter cf> as 
a function of fi/U obtained by the perturbative MF method. 
We set zt/U — 0.13. The thick solid, thin solid, and dashed 
lines represent the ground, metastable, and unstable states, 
respectively. In the regions (I- IV), the energy function has 
the shape shown in Figs. |15f I-IV) , respectively, (b) The cor- 
responding hysteresis loop structure when cycling the value 
of n/U. 



SF solution) in the intermediate regions (II) and (III). 

The ground state of regions (I) and (II) is a SF state, 
which is located at the global minima of E((f>). If one 
increases fi/U from the initial SF state, the system re- 
mains in the metastable SF phase even in region (III) 
since E{(j>) still has local minima at finite <j>. The SF 
state is destabilized in region (IV) and the transition to 
the MI phase occurs due to the disappearance of the lo- 
cal minima [sec the illustrations in Fig. [TTlb)]. On the 
other hand, when starting from an initial MI state in re- 
gion (IV), it remains locally stable until region (II) and 
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FIG. 18: (a) The same as in Fig. UTTa) for zt/U = 0.153. 
(b) The anomalous behavior in the hysteresis process. The 
system keeps in the SF phase ignoring the solutions of MI 
states. 



then the system undergoes a transition to the SF phase 
in region (I). As a result, the hysteresis cycle forms a loop 
structure as shown in Fig. [T7TbV 

For zt/U close to the tip of the MI lobe, the hysteresis 
process exhibits the anomalous behavior as discussed in 
Sec. IIII Al The upper panel of Fig. [T51a) shows the solu- 
tion curve of as a function of \x/U for zt/U = 0.153. 
The full MF result shown in Fig. [TTTa) can be repro- 
duced almost perfectly by the perturbative MF calcu- 
lation with up to the sixth-order perturbations. In this 
case, the system can undergo the transition only unidirec- 
tionally from the MI to SF phase as shown in Figs. lTlT b) 
and [TITc) . The cause of this anomalous hysteresis can 
be understood from the behaviors of the coefficients in 
the Ginzburg-Landau energy plotted in the lower panels 
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of Fig. [IBTa). Unlike Fig. UTTa). the curve of a\ — Aa^a^ 
crosses twice since the energies of SF and MI become 
identical at the two points. The curve of a-i also takes 
twice at the points which correspond to the metasta- 
bility limits of the MI phase. However, the curve of 
a\ — 3a2a6 does NOT touch 0, which means that the 
coefficients always satisfy the inequality a^ < —yJ3a2a,Q 
in the anomalous hysteresis. This fact is essential for the 
emergence of the anomalous hysteresis. The behaviors of 
the coefficients indicate that the energy profile changes 
as (I) — ^(11) — ^-(IH) — ^(11) — ^(1) without going through re- 
gion (IV) when [ijU increases. Therefore, E{4>) always 
has minima at </> ^ and a stable SF state exists over 
the entire range of /J./U. In consequence, the system is 
caught in a local minimum of E(4>) at </> ^ and the 
SF state is not destabilized dynamically [sec the illustra- 
tions in 118( b)]. This is the cause of the unidirectional 
transition process in the anomalous hysteresis. 

When Ui y^ U2, the sixth-order Ginzburg-Landau en- 
ergy is given as a function of (f>\ and (f>2- The energy 
function E((f>i,(j)2) consists of many terms proportional 
to oc <j)\ m <j)l n (0 < m + n < 3). When \Ui - U 2 \ < Ui 
(e.g., in the case of Scc. lIHB]) , the values of <\>\ and </>2 are 
almost identical and the first-order SF-to-MI transitions 
should be also described by the sixth-order Ginzburg- 
Landau energy in a similar way. 
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zt/U 

FIG. 19: The ground-state phase diagram of the spin-1 Bose- 
Hubbard model within the MF approximation for U2/U0 — 
0.04. Second- and first-order phase transitions are indicated 
by the thin and thick lines, respectively. The dashed (dash- 
dotted) lines represent the metastability limits of SF (MI) 
phase, and the dots mark the tricritical points, where the 
transition changes from first to second order. The inset is an 
enlarged view of the region around the tip of the p = 2 MI 
lobe. 



IV. SPIN-1 BOSONS 



In Sec. IIIID1 we have shown that the anomalous 
hysteresis can be qualitatively understood from the 
Gintzburg-Landau theory, in which the dependence of 
the coefficients on microscopic system parameters plays 
an essential role in causing the unidirectional behavior. 
The anomalous hysteresis is now mathematically well es- 
tablished as another type of first-order transition phe- 
nomena different from the conventional hysteresis-loop 
formation. However, such an anomalous hysteresis has 
never been observed in experiments on either cold-atom 
systems or solid-state materials. Therefore, in this sec- 
tion we add another example to the list of possible sys- 
tems that exhibits the anomalous hysteresis in order to 
provide a further guideline for achieving the first obser- 
vation in experiments. 

Here we consider the system of spin-1 bosons confined 
to an optical lattice [Mill [HH86J]- The previous stud- 
ies [38l442l ] have predicted that the phase transition be- 
tween the SF and MI states can be of first order at even 
fillings. We model the spin-1 bosons in a hypercubic 
lattice with the following spin-1 Bosc- Hubbard Hamilto- 



H = -tJ2 E ®,A* +H.c)-/i^i 



(j,l) o-=0,±l 

3 3 



*-2n j ),{22) 



b- a creates a boson with spin a at site j 



of 



where b] 

a D-dimensional hypercubic lattice, hj = J2a- 

?3 = E^LF-A*': and ¥ aa , = (F^Fv',FZ a> ) 

I I I I 

are the standard 3x3 spin-one matrices [40J, [83J . This 
system can be realized using a gas of alkali atoms such as 
23 Na, 39 K, and 87 Rb with hyperfine spin F = 1 [87rll. 
The strength of the interactions Uq and U% are related 
to the scattering lengths cio and 02 for 5 = and 5 = 2 
channels, and it is known that 23 Na has antiferromag- 
netic interaction U2 > while 87 Rb is ferromagnetic, 
U2 < 0. For antifcrromagnctic interaction Ui > 0, the 
SF-to-MI transition can be of first order near the tips of 
the MI lobes at even filling factors [H El ■ 

Figure Q1J] is the ground-state phase diagram of the 
spin-1 Bose-Hubbard model for U2/U0 = 0.04 obtained 
by the single-site MF approximation with the decoupling 
similar to Eq. ([5]). This value of U2/U0 corresponds to the 
experiment with 23 Na atoms 89M9ll|. For spin-1 bosons, 
one has to solve the self-consistent equations for three 
mean fields 4> a = (b a ) (a = 0, ±). It is found that the SF 
region has a polar solution 0o 7^ 0, 4>± — as the ground 
state [83J. The phase boundaries between SF and MI 
obtained by the MF decoupling procedure are identical 
to the ones obtained before by the Gutzwiller approxi- 
mation [401 due to the equivalence of the two single-site 
approximations. In addition to the second- and first- 
order boundaries, we also plotted the metastability limits 
of SF and MI, at which the corresponding local minima 
of the energy function disappear. When the filling fac- 
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FIG. 20: Finite-temperature phase diagram of the spin-1 
Bose-Hubbard model with U2/U0 = 0.04. The chemical po- 
tential p is tuned to be the value at the tip of the p = 2 MI 
lobe; p/Uo = 1.338. Second- and first-order phase transitions 
are indicated by the thin and thick lines, respectively. We also 
plot the value of the Mott gap Ami/£/o- Lower panel: the en- 
larged view of the upper panel with adding the metastability 
limits the SF and MI phases. 



tor p = (n) is even, the MI state consists of p/2 singlet 
pairs with S = while it has spin S = 1 at odd fill- 
ing factors 



Only the case of even filling factors can 
exhibit first-order SF-to-MI transitions. One has to go 
beyond the simple single-site approximation in order to 
discuss more detailed magnetic structures inside the MI 
phases dHH. 



The phase diagram with the metastability limits has a 
very similar structure to the one of the two-component 
Bose-Hubbard model with repulsive inter-component in- 
teraction in Fig. [8] As shown in the finite-temperature 
phase diagram (Fig. I2TJ)) . the transition between the SF 
and p = 2 MI phases remains first order until T/U\ rj 
0.11 and there is no quantum critical point. Therefore, 
the system of spin-1 bosons in an optical lattice is an- 
other qualified candidate for simulating first-order transi- 
tion phenomena in cold-atom systems, although it needs 
to be cooled to sufficiently low temperatures. We can 
see that the phase diagram in Fig. [8] has the geome- 
try required for the anomalous hysteresis to appear; the 
Mott lobes at even fillings are surrounded by first-order 
transition boundary with the SF phase. Therefore, the 
unidirectional behavior in the hysteresis is expected to 
be observed also in the spin-1 system upon varying the 
chemical potential h/Uq. To obtain the anomalous hys- 
teresis, the optical-lattice depth has to be tuned so that 
the value of zt/Uo lies within <5 a shown in the inset of 
Fig. EH 



V. SUMMARY 

In summary, we have studied first-order transition phe- 
nomena of Bose gases trapped in optical lattices, espe- 
cially focusing on an anomalous hysteresis behavior that 
was proposed in Ref. |52J. We first reconsider the ground- 
state phase diagram of the hardcore bosons with full 
dipolc-dipole interaction in a triangular lattice [5l|, [52[ 
by means of the cluster mean-field method with cluster- 
size scaling [571 ] . Our scaled phase diagram is in good 
agreement with the numerical deta obtained by quantum 
Monte Carlo calculations [5l| within the error bars, ex- 
cept that our result predicts the existence of several addi- 
tional solid phases. We demonstrated that the first-order 
melting transition can exhibit an anomalous hysteresis 
behavior upon sweeping the value of the chemical poten- 
tial, as in the truncated model with only nerest-neighbor 
interaction [521 ] . In the anomalous hysteresis cycle, once 
a solid state melts into superfluid, it cannot be solidified 
back again. Unlike in a conventional first-order transi- 
tion, the hysteresis curve does not form a hysteresis-loop 
structure. We have also studied binary Bose mixtures 
with inter-component repulsion loaded into a hypercu- 
bic lattice. We found a similar unidirectional behavior 
in the superfluid- Mott insulator transition of this system 
upon sweeping the value of the chemical potential; the 
transition from superfluid to Mott insulator cannot oc- 
cur whereas the transition in the opposite direction is 
allowed. 

A common feature of the above two bosonic systems 
is a chractcristic geometry of the phase diagram; in both 
systems, a phase region of lobe shape (the solid phase for 
the former case and the Mott insulator for the latter case) 
is surrounded by the first-order boundary. This geome- 
try is indeed the necessary condition for the anomalous 
hysteresis to emerge. As another example that has such a 
geometry, we also addressed the system of spin-1 bosons 
in a hypercubic optical lattice. Moreover, in our recent 
work [96f . we found that a spin-dimer system on a tri- 
angular lattice can also exhibit the anomalous hysteresis 
behavior in the first-order magnetic transition induced 
by controlling the magnetic field instead of the chemical 
potential. Having these examples, we conclude that the 
anomalous hysteresis is a ubiquitous phenomenon of sys- 
tems with a phase region of lobe shape that is surrounded 
by the first-order boundary. 

Furthermore, taking the system of binary Bose mix- 
tures with equal intra-componcnt interactions as a simple 
example, we constructed a sixth-order Ginzburg-Landau 
theory for the anomalous hysteresis in order to mathe- 
matically establish the unconventional first-order tran- 
sition. Using the perturbative mean-field method, we 
determined the coefficients 02, 04, and ag as functions of 
the microscopic parameters t, U, C/12, and pL. We found 
out that in the anomalous hysteresis, a^ never exceeds 
— V / 3o206 an d thereby the energy landscape always has 
minima at finite superfluid order parameter. Therefore, 
superfluid states remain stable over the entire hysteresis 
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region, resulting in the anomalous unidirectional hystere- 
sis. 

Assuming binary mixtures of 87 Rb atoms in the two 
hyperfine states [3HZ| confined to a simple-cubic optical 
lattice, we estimated the required range of the tuning 
of lattice depth Vq for a clear observation of first-order 
transition phenomena. It is expected that the estima- 
tion shown here should stimulate further experiments in 
this direction. We stress that the observation of the first- 
order transition will significantly expand the applicabil- 



ity of ultracold atomic/molecular systems as a quantum 
simulator of strongly interacting many-body physics. 
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